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On s'interesse ici a des problemes ou I'inconnue est une fonction : 
9 Chute libre : Trouver la vitesse v(t) satisfaisant 

mv'(t) = mg- kv{tf Vf 

9 Pendule : Trouver Tangle 9{t) satisfaisant 

ml 9"(t) = -mg sin{9(t)) - lc f 6'{t) Vf 

9 Masse-ressort : Trouver le deplacement x(t) satisfaisant 

x"(t) + cx'(t) + oj 2 x{t) = Vt 

9 Croissance de population : Trouver le nombre d'individus N(t) 

N'(t) = (n^-m^N(t))N(t) Vf 
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Introduction Forme generale 



Equation differentielle 



Une equation differentielle est une relation faisant intervenir les 
derivees d'une fonction inconnue y(t) a determiner sur I'intervalle 
[a,b]: 

F(t,y(t),y'(t),y"(t),...,y^(t)) = Vte[a,b] 



Definition 1 



L'ordre d'une equation differentielle est determine par d'ordre 
maximale de la derivee de la fonction inconnue (variable dependante). 



Exemple 



a Ordre 1 : y'(t) = ty{t) alors F{t,y,y') = y'(t) - ty(t) 
9 Ordre 2 : y"{x) = y'(x) - y(x) alors 

F(x,y,y',y") = y"(x)-y'(x) + y(x) 
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9 Contrairement a une equation algebrique, I'inconnue du probleme 
est une fonction. 

» La solution generale depend de un ou plusieurs parametres. 
9 Pour avoir unicite, il faut preciser des conditions 

« initiales (valeurs connues pour la plus petite valeur de la variable 

independante), generalement la variable independante est associe 

au temps (et notee f) 

y'(t) = ty(t) y(0) = 1 

» des conditions aux limites, (valeurs connues a une/plusieurs 
extremites de I'intervalle de la variable independante), 
generalement la variable independante est associe a une variable 
d'espace (et habituellement note x) : 

y"(x) = 2a y(0) = c, y(1) = a + £> + c 
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Equations differenti 



Equation d'ordre 1 



La formulation generale d'une equation differentielle d'ordre 1 avec 
condition initiale est : 



y'(t) = f(t,y(t)) 

y{k) = yo 



Exemple 1-2 



• y'(t) = t avec y(0) = 1 donnant y(t) = f + 1 

( 2 -3 

9 y'(t) = ty{t) avec y(3) = 42 donnant y{t) = 42.e~ 



9 En general il est impossible de resoudre analytiquement ces 
equations. 

9 Dans les cas simples on peut faire une analyse des directions de 
I'equation donnant une "image" de la solution 
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Introduction Equations differentielles d'ordre 1 



Exemple : Champs des directions de I'equation differentielle 



At) 

y(o) 



-y + t + ^ 



1 



Dans le plan (t, y) on fait une 
analyse du "champ vectoriel" : 
Pour chaque point (t, y) on 
peut calculer la derivee de y. 
Exemple en (1,1) on a 
y'(1) = 1 eten (-1,1) on a 

On peut alors "remplir" le plan 
et voir "apparaTtre" une solution 



\\\ \V\iVV 

\ \ \ V\\\V- 
\ \\ Si \ \n- 

s,w^.-« — - 
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Equations different 



Rappels 



9 On s'interesse aux methodes numeriques de resolution. 

9 Au niveau numerique, on ne peut pas evaluer la solution y(t) en 
tout point t. 

9 On va chercher une apprioximation de la solution pour certaines 
valeurs de la variable independante t qu Ton notera t,. 



Vocabulaire 



9 La distance entre deux valeurs successives f,- +1 et l, est note 
hj = tj + -\ - tj s'appelle le pas de temps 

9 On utilisera generalement un pas de temps constant note h 

9 On note y(f,-) la solution analytique evaluee en t = t-, 

9 On note y ; la solution approximative en t = t-, obtenue par une 
methode numerique 
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Introduction Discretisation 



Discretisation 



Supposons que Ton desire calculer la solution y(t) d'une equation 
diffrentielle pour t > t. Une discretisation de I'equation definissant 
y{t) est definie comme 

• la donnee d'un ensemble discret de valeur pour la variable 
independante : en general, on fixe le pas h et on pose 



ti = t + ih V/' = 0,1. 



N 



Produisant la suite 

t < U < t 2 < ■ ■ ■ < ti < ■ ■ ■ < Nh 

9 et d'une methode numerique (schema) fournissant une 
approximation de la solution exacte au points t-, : 

y, « y{tj) 
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Soit le probleme 

f y'(t) = f(t,y(t)) 

I y(k) = yo 

On se donne un pas de temps h et t; = t Q + ih. 

Pour approximer /(^ ) = y(t Q + ft) on propose d'utiliser I'information 
connue : 

to,y(t ) ety'(fo) = f(f ,y(f )) = f(t ,y ) 

On construit la droite d Q (t) passant par (t ,y ) et de pente f(t ,y ) : 

d (0 = yo + (f-Wfo,yo) 

et on propose y : = d {h) = y + {h - t )f{t ,y Q ) = y + hf{t ,y ) 
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Pour approximer y(fe) = y(k + 2h) on utilise le meme procede. On 
construit une droite avec I'information connue en t\ : la droite d^(t) 
passant par (^ , yi ) et de pente f(U , y-\ ) : 

tfi(0 = yi+C- Wi,yi) 

et on propose y 2 = d, (fe) = yi + (fe - fi Wi , /i ) = /i + W(fi , yi ) 
et on recommence au pas suivant donnant 

yi+A=yi + hf(t h yi) i = 0,... 



Definition 



On qualifie d'explicite un schema ou revaluation de I'approximation 
ne necessite pas de resolution d'un systeme non-lineaire. Plus 
simplement I'approximation s'exprime en fonction de valeurs connues 
uniquement. 
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On a 



y'(t) = f(t,y(t)) Vt€[t ,Nh] 



y , (t i ) = f(t i ,y(t i )) V/ 



On arrive a la meme methode en considerant I'approximation de la 
derivee y'(tj) par la formule de difference avant 



W,y/) «/('/) 



y/+i - y/ _ y/+i - y/ 



menant a 



y+i - y 



W,//) v/ 



Cette methode est connue comme la methode de Euler explicite. 
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Methode d'Euler explicite 



• Etant donne un pas de temps h, une condition initiale (fo,yo) et un 
nombre maximal d'iterations N 

9 i = 

9 y/+i =yi + hf(t,,yi) 
9 f /+1 =ti + h 
9 Ecrire f, +1 et y, +1 
9 Si / = A/ arret 
Q / = / + 1 retour 1 
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Erreur 



II faut cependant noter qu'en general 

y, ± y(ti) et y'(t,) = f{t h y(t,)) + f(t,,y,) 



done 



di(t)^y(ti) + (t-t,)f(ti,y(ti)) 



Mais alors I'erreur commise dans I'approximation de y, est introduite 
dans I'approximation de y /+1 et on a propagation de I'erreur ! 
La propagation de I'erreur n'est pas exclusive a la methode de Euler 
explicite. Elle est typique des methodes numeriques pour les eq. diff. 
On aura propagation de I'erreur en passant d'une iteration a I'autre, et 
de maniere generale I'erreur d'approximation 

\y(ti)-yi\ 

augmente (legerement, quand tout va bien) avec /'. 
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Methode d'Euler explicite Methode d'Euler : exemple 



Soit I'equation differentielle : 

y'(t) = -ytf + t + v, 
y(0) = 1. 

qui admet la solution exacte y(t) = e~ f + t. 



Le schema d'Euler s'ecrit 

yo = 

y n +i = 

Pour h = 0.1 on obtient 

y =y + 0.1(-y + ? + 1) = 
y 2 = yi +0.1(-y +t^+^) = 
y 3 = y 2 + 0.1 (-y 2 + 0.2 + 1) 



1, 

yn + h{-y n + t n + }). 



1 +0.1(-1 +0 + 1) = 1 

1 + 0.1 (— 1 +0.1 +1) = 1.01 

= 1.01 +0.1(— 1.01 +0.2 + 1) 



1.029 
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Methode d'Euler explicite Methode d'Euler : exemple 



Methode d'Euler :y'(t) = 


y{t) + t + i 


tn 


Y(tn) 


Yn 


Wn)-y n \ 


0.0 


1.000000 


1.000 000 


0.000000 


0.1 


1.004 837 


1.000 000 


0.004 837 


0.2 


1.018731 


1.010 000 


0.008731 


0.3 


1.040818 


1.029 000 


0.011818 


0.4 


1.070302 


1.056100 


0.014220 


0.5 


1.106531 


1.090490 


0.016041 


0.6 


1.148812 


1.131441 


0.017371 


0.7 


1.196585 


1.178 297 


0.018288 


0.8 


1.249329 


1.230467 


0.018862 


0.9 


1.306570 


1.287420 


0.019150 


1.0 


1.367 879 


1.348 678 


0.019201 
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Quel est I'effet de h sur I'erreur \y(t n ) - y n \ ? 



Definition 2 



Une methode de resolution est dite methode a un pas si elle est de la 
forme 

y n +i =yn + h<f>(t„,y„) 

ou cj) est une fonction quelconque. La methode est a un pas si elle 
exige uniquement la solution au temps precedent pour la resolution au 
temps courant. 

Une methode est a pas multiples si la solution au temps t n+ -\ exige la 
solution a plusieurs valeurs de temps precedent t n , ? n -i . tn-2,-- 
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Definition 3 



Erreur de troncature locale au temps t n est 

y(f„+i ) - y(t„) 



'n+1 



(h) 



%,y(t n )) 



Remarque 



On utilise la solution exacte dans I'expression de I'erreur de troncature 
locale. C'est parce qu'on veut mesurer I'erreur introduite par 
I'utilisation du schema numerique pour un pas donne en supposant la 
solution exacte aux iterations precedentes. 

y(f n+1 ) - y n+ i = y(f„+i )-y„- h<j>(t n , y n ) 
En supposant la solution exacte en t n on a 

y(f„+i ) - y n+ i = y{t n +i ) - y{t„) - h<f>{t„,y(t„)) = /7r n+1 (h) 
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Revenons a la methode d'Euler explcite. Dans ce cas <p(t,y) = f{t,y). 
En utilisant Taylor (encore) autour de t n on a 

y(f n+1 ) = y(t n + h) = y{t n ) + y'(t n )h + t!Mt + (h 3 ) 
= y(t n ) + f(t n , y(t n ))h + ^IMt + (h 3 ) 



alors 



y(tn+i) - y{t„) 



f(t„,y{t„)) + 



2 

y"(t n )h 



+ 0(/7 2 ) 



L'erreur de troncature locale est done : 

r fM y(t n+ i) - y(t n ) .. . „ y"{t n )h 2 



C'est-a-dire 



r n+1 (/7) = 0(/7) 
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Conclusions 

9 Si on suppose que y n correspond a la solution alors sur chaque 
intervalle [t n , t n+ -\], I'erreur E n = \y(t n ) - y n \ sera de I'ordre de 

0(/7 2 ). 

9 Cependant par la propagation des erreurs on a au mieux une 
erreur correspondant a I'erreur de troncature locale, i.e. une 
erreurd'ordre 0(h) 

Remarques 

9 L'ordre de la troncature locale donne I'ordre de la methode 
numerique, la methode de Euler explicite est d'ordre 1 . Attention 
de ne pas confondre avec I'ordre de I'equation differentielle. 

9 La methode de Euler explicite est facile d'emploi mais peu 
precise. C'est pourquoi elle est peu utilisee. 

9 Pour reduire de moitie I'erreur dans la methode de Euler explicite 
on doit diviser le pas de temps par 2 puisque la methode est 
d'ordre 1. 
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Methode d'Euler explicite Exemple, erreur 



Revenons a I'exemple 3 : 

y'{t) = -y{t) + t+A y(0) = 1 y(t) = e- t + t 

En utilisant Euler explicite avec h = 0.1 et 0.05 on trouve 

tj /? = 0.1 erreur h = 0.05 erreur 

0.1 1.0 0.004837 1.007375 0.002537 

0.2 1.01 0.008731 1.0237809 0.005050 

On a bien un reduction de moitie de I'erreur pour une reduction de 
mopitie du pas de temps. 

Pour gagner en precision on peut "etendre" la methode de Euler 
explicite en utilisant le developpementde Taylor. Cette approche 
permet d'avoir une erreur de troncature locale d'ordre plus elevee. 
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On utilise le developpement de Taylor pour la fonction y(t) 

y"{t n )rf 



y{t n+ , ) = y(t n + h) = y(t n ) + y'(t n )h + 



+ 0(/? 3 ) 



Or 



y"(t) 



df(t,y(t)) , 8f(t,y(t)) „, 



at 



+ 



dy 



y'(t) 



df(t,y(t)) , 8f(t,y(t)). 



On obtient 

y(t n+ ,) = y(t n ) + hf(t n ,y(t n )) 
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En negligeant les termes d'ordre superieurs a 2 en h, et en remarquant 
que le membre de droite est compose d'elements connus (on peut 
evaluer les derivees de f et elles sont evaluees a des points connus). 

On obtient : 



y(t n+ :)^y(tn) + hf(t n ,y(t n )) 

h 2 fdf{t n ,y{tn)) 



+ 



ot 



+ 



df(t n ,y(tn)) 

dy 



f(tn,y(tn)) 
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On "remplace" la solution exacte y(t n ) par I'approximation y n et on a le 
schema : 

Xn+1 - /n + nT{tn,yn) + ~2 I ^ 1 ^ i[}n,yn) 

Puisqu'on utilise y n plutot que la solution exacte on aura encore 
propagation d'erreur d'un pas de temps a I'autre. 



Erreur de troncature locale 



Suivant la notation on a 

Part construction (on peut facilement montrer que) 

r n+1 (/7) = 0(/7 2 ) 

La methode est d'ordre 2. 
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Schema de Taylor d'ordre 2 J 


Q Etant donne un pas de temps h, une 


condition initiale 


(k 


yo) 


et un 


nombre maximal d'iterations N 










Q PourO < n < N : 










* ,-<, x h 2 fdf(t n ,y r 
© y„+i = y n + hf(t n ,y n ) + y f — ^ 


) , df(t n ,y n ) fft 
1 dy t{tr 


,y r 


') 




O ? n+ 1 =t n + h 










Q Ecrire f n+1 et y n+1 










Q Arret 
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Methode de Taylor : y'{t) = -y(t) + t+*\ 


ti 


y(ti) 


// 


\y(ti)-yi\ 


0.0 


1.000000 


1.000 000 


0.000000(0.000000) 


0.1 


1.004 837 


1.005 000 


0.000163(0.004 837) 


0.2 


1.018731 


1.019 025 


0.000294(0.008731) 


0.3 


1.040818 


1.041218 


0.000400(0.011818) 


0.4 


1.070302 


1.070 802 


0.000482(0.014220) 


0.5 


1.106531 


1.107 075 


0.000544(0.016041) 


0.6 


1.148812 


1.149404 


0.000592(0.017 371) 


0.7 


1.196585 


1.197210 


0.000625(0.018288) 


0.8 


1.249329 


1.249 975 


0.000646(0.018862) 


0.9 


1.306570 


1.307 228 


0.000658(0.019150) 


1.0 


1.367 879 


1.368 541 


0.000662(0.019201) 
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9 equation differentielle : unicite de la solution obtenue en ajoutant 
des conditions initiales ou au bords i.e. on se donne la solution en 
certains points. 

9 I'ordre de lequation : determinee par le degre de la plus haute 
derivee de I'inconnue 

9 le pas de temps, discretisation et schema numerique, schema 
explicite 

9 schema d'Euler explicite : 

y n +i =yn + hf{t n ,y n ) 
9 Dans les schemas pour les eq. diff. I'erreur se propage 
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9 Schema a un pas general : 

y n+ i =yn + h(f>{tn,y n ) 

9 Erreur de troncature locale pour un schema qui donne I'ordre du 
schema. 

9 Pour augmente I'ordre de I'erreur de troncature on produit la 
methode de Taylor d'ordre 2 : vient du developpmenet de Taylor 
de y autour du point t n avec troncature des termes d'ordre 3 en h. 

v -v+hf(t y\ I ^ f df (tn,yn) df(t n ,y n ) 

yn+1 - y n + nt[t n ,y n ) + — I — 1 — T(t n ,yn) 

9 Euler explicite est un schema d'ordre 1 
9 La methode de Taylor est d'ordre 2. 
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Revenons a y'(t) = -y(t) + 1 + 1 . En utilisant Euler explicite avec 
h = 0.1 on a 

|y(1.0)-y 10 | =0.019201 

et avec la methode de Taylor on a 

1/(1-0) -/io| =0.000662 

Avec le schema d'Euler, quelle valeur de h faut-il pour avoir une erreur 
en t — 1.0 comparable a I'erreur due a la methode de Taylor? 

Puisque Euler est un schema d'ordre 1 on a 

|y(1.0)-y n |«C/7 

On veut une erreur qui soit 0.01 9201 /0. 000662 « 30 plus petite alors 
il faut prendre h a peu pres 30 fois plus petit ! 
On devra done prendre h m 0.1/30 
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Conclusions 

9 En prenant h rj 0.00333 on aura une erreur comparable 

CEPENDANT on devra faire 300 iterations pour atteindre la valeur 
de t = 1.0 au lieu de 10. 

9 On peut done controler I'erreur jusqu'a un certain point en 
modifiant le pas de temps mais e'est au detriment du nombre 
d'iterations. 

9 Le tableau p 355, procede au meme traitement mais avec la 
methode de Taylor d'ordre 2. En faisant une division par 2 du pas 
de temps on reduit par un facteur 4 I'erreur . On note qu'en fait le 
facteur 4 s'obtient quand h est tres petit. 

9 Tout comme pour le schema d'Euler le nombre d'itetations grandi 
(multiplie par 2 a chaque fois). 
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Le choix du pas de temps semble arbitraire. Qu'en est-il ? 
Considerons le probleme 



At) 

y(o) 



-5y(f) 



1 



Dont la solution est y(t) = e~ 5t . En utilisant Euler explicite avec un pas 
de temps h donne on a 

/o = 1 y n +i =y n - 5hy n = y„(1 - 5/7) n > 

Mais alors 

y n = y (1 - 5/7) n = (1 - 5/7) n 

II est claire que dans le cas ou (1 - 5/7) < la solution oscillera entre 
des valeurs positives et negatives : 

/7 = 0.25^ y n = (-0.25)"^ 
y = 1,/1 = -0.25, y 2 = 0.0625, y 3 = -0.015625... 
Mais on sait deja que la solution est strictement positive ! 
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Conclusion 

9 Le schema produira une approximation valable pour certaine 
valeur du pas de temps seulement ! 

9 Dans le cas present on veut h < 0.2. 

9 Puisqu'en partie la construction des schemas ainsi que I'erreur de 
troncature qui "qualifie" nos schemas est basee sur le 
developpementde Taylor, la conclusion habituelle s'applique. II 
est dangeureux d'utiliser des points trop "loin" les uns des autres : 
il faut des pas de temps petits si on veut etre prudent ! 
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Puisque modifier le pas de temps n'est pas ideal pour ameliorer la 
qualite de la solution, on regarde ce qui se passe en augmentant 
I'ordre de la methode de Taylor (en tronquant "plus loin") : 

Par exemple pour avoir un schema d'ordre 3 : 
y(t n+ ,)=y(t n ) + hf(t n ,y(tn)) 

h 2 fdf(t n ,y(t n )) , df(t n ,y(tn)) 



+ y 



+ 



+ 
+ 
+ 



dt dy 

/7 3 f d 2 f(t n ,y(t n )) d 2 f(t n ,y(t n )) 
6 V dt 2 dtdy 

df(tn,y(tn))df(t n ,y(tn)) , d 2 f(t„, y(t n) ) 

dt dy 

df(tn,y(tn))df(t n ,y(tn)) 



+ 



f(tn,y(tn)) 
f(tn,y(tn)) 

(f(tn,y(tn))Y 



dy 



dy 



dy 2 
(f(t n ,y(t n ))) 2 ) fO(A7 4 ) 



Donne un gain en precision au detriment de la simplicite de 
I'expression a evaluer. 
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On voudrait une/des methodes d'ordres de plus en plus grand tout en 
evitant la complexite devaluations des methodes de Taylor d'ordre 
elevee. 

Dans les methodes de Taylor on augmente I'ordre en "extrayant de 
I'information" des derivees de la fonction f (c-a-d y'). C'est aussi la 
source des difficultes car on devra "complexifier" le schema en 
evaluant des derivee d'ordre de plus en plus eleve. 

Pour avoir des schemas simples on evitera d'utiliser les derivees de la 
fonction /. Pour garder I'ordre des methode de Taylor sans utiliser les 
derivees il nous faudra "prendre I'information ailleur" mais il ne reste 
plus que f ! On devra done faire des evaluations de f en plusieurs 
points. 

C'est la premisse des methodes de Runge-Kutta (RK). 
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Methodes de Runge-Kutta Methode de Runge-Kutta d'ordre 2 



On s'interesse d'abord a I'ordre 2. 

La base de la methode de Taylor d'ordre 2 est le developpement : 

y(t n+ :)=y(tn) + hf(tn,y(tn)) 

+ *(%&&» + ^0m f{tnMtn) ^ + 0(/?3) 

On veut un schema de meme ordre de precision et ne faisant pas 
apparaTtre de derivees. On propose comme base des schemas d'ordre 
2 : 

yU n +i ) = y(t n ) + ai hf(t n , y{t n )) + a 2 hf{t n + a 3 h, y(t n ) + a 4 h) 
Ou a-\ , a 2 , a 3 et a 4 sont a determiner (pour satisfaire I'ordre 2). 
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Pour trouver ces coefficients on fait le developpement de Taylor (cette 
fois en 2 variables) de la fonction f : 

f{t n + a 3 h,y{t n ) + a 4 h) = f(t n + v ,y{t n ) + v^ ) 

_ , ( t,, y(W) + Vo ?i^mi + ^mgm + 0( ^ 



dt 

df(t n ,y(t n )) 



dy 

df(t n ,y(tn)) 



f(t n , y(t n )) + a 3 h -^'" d y'"" + a A h -^" d '^"> + 0(h 2 ) 

On remplace ensuite dans I'expression : 
/(fn+i) = y(t n ) + {a^ + a 2 )hf{t n ,y(tn)) 



9f 



<9y 



II ne nous reste plus qu'a identifier les termes en les comparant au 
termes utilises pour la methodes de Taylor. 
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Methodes de Runge-Kutta Methode de Runge-Kutta d'ordre 2 



y(t n+ :)=y(tn) + hf(t rh y(t n )) 

h 2 fdf(t„,y(t n )) , df(t n ,y(tn)) 



+ y 



at 



+ 



dy 



f(t„,y(tn)))+0{h a ) 



y{t n+ \) = y(t n ) + (ai + a 2 )hf(t n ,y(t n )) 

+ £ ( 2 a 2 a 3 ^fM + 2a 2 a 4 «M ) + O(tf) 

Pour que les deux expression concordent il faut : 

a! + a 2 = 1 

1 



32^4 



32^3 — p 

'(fti,y(f n )) 
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Methodes de Runge-Kutta Methode de Runge-Kutta d'ordre 2 



y{t n +i ) = y(tn) + ai hf(t n , y{t n )) + a 2 hf{t n + a 3 h, y(t„) + a 4 h) 

1 f(tn,y(tn)) 



a-\ + &2 — 1 , 32^3 



32^4 



2' 2 

Le systeme a resoudre est sous-determine (plus d'inconnues que de 
variables) on est done assure d'une solution et permet differentes 
variantes. 

Deux variantes "classiques" menant a des schemas frequemment 
rencontre : 

9 a 1 = a 2 = \, a 3 = 1 et a 4 = f(t n ,y(t n )) mene au schema d'Euler 
modifiee en remplagant y(t n ) par y n et en negligeant les termes 
en 0(h 3 ) dans I'expression 

• a^ = 0, a 2 = 1 , a 3 = \ et a 4 = ^ f "'^ fn)) mene a la methode du 
point milieu en remplagant y(t n ) par y n et en negligeant les 
termes en 0(h 3 ) dans I'expression 
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Schema RK2 : Euler modifiee 


9 


On se donne h, {t Q ,y Q ) et un nombre maximal d'iterations N 


9 


Pour < n < N : 




o 


y = y n + a 4 h = y n + hf(t n , y n ) 




o 


y n+1 =y n + a : hf(t n , y n ) + a 2 hf{t n + a 3 h, y n 

= yn + ^(f(tn,yn) + f(tn + h,y)) 


+ a 4 h) 


e 


t n+ -\ =t n + h 




O Ecrire f n+1 et y n+1 




Q 


n = n+1 retoura 1 





L'etape 1 correspond a la mehode de Euler explicite, d'ou le nom de 
methode d'Euler modifiee. La methode est en fait composee de deux 
calculs : une prediction par l'etape 1 suivie d'une correction 
(amelioration) de cette prediction dans l'etape 2. C'est une methode 
du type predicteur-correcteur 
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Methodes de Runge-Kutta Methode du point milieu 



Schema RK2 : Methode du point milieu 



• Etant donne un pas de temps h, une condition initiale (fo,yo) et un 
nombre maximal d'iterations N 

9 PourO < n < N : 

O /ci =a 4 h = hf(t n ,y n ) 

9 y„+i = y n + a : hf(t n , y n ) + a 2 hf{t n + a 3 h, y n + a 4 h) 



yn + hf(t n + ly n + ^ 



O t n+: =t n + h 

O Ecrire t n+ : et y n+1 

Q n = n+ 1 retoura 1 



On evalue f au point t n + 5 d'ou le nom de la methode. 
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Methodes de Runge-Kutta RK d'ordre > 2 



Le procede de construction utilise pour I'ordre 2 peut etre applique 
pour les ordres suivant. II faut noter que cette approche devient de 
plus en plus complexe. 

Dans tous les cas on obtient un systeme sous determine a resoudre 
pour les coefficients. 

Le schema vedette de la famille des RK est le schema d'ordre 4 
suivant : 

1 / h 

y n +i =Yn+Q (hf(t n ,y n )+2hf{tn+^,yn+hf{tn,y n ) 

hi hi 

+ 2hf{t n + ^y n + hf(t n + ^,y n + hf{t n ,y„))) 

hi hi 

+hf{t n + h,y n + f(t„ 4 -,y„ + hf{t„ +^,y n + hf{t„,y„)))) 

Ce schema est tellement utilise qu'il est connu sous le nom de "RK4" 
oubliant qu'en fait il y a une "famille" de methode de Runge-Kutta 
d'ordre 4. = i ^<^ 
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m 



Schema RK4 ^ 


O Etant donne un pas de temps h, une condition 


initiale (fo>yo) et un 


nombre maximal d'iterations N 




Q PourO < n < N : 




O /ci =hf(t n ,y n ) 




Ok 2 = hf[t n + ^y n + ^ 




Q k 3 = hf(t n + ~,yn + 'j) 




O k 4 = hf{t n + h,y n + k 3 ) 




© y n +i = y n + g (/ci + 2/c 2 + 2/c 3 + /c 4 ) 




O t n+: =t n + h 




O Arret 





Maths-Info (ENSAIV 



Methodes de Runge-Kutta Methode de Runge-Kutta d'ordre 4 



Methode de Runge-Kutta : y'(t) = -y(t) + t+*\ 


ti 


y(ti) 


// 


\y(ti)-y,\ 


0.0 


1.0 


1.0 


0.0 


0.1 


1.004 837 4180 


1.004 837 5000 


0.819 x 10" 7 


0.2 


1.018 730 7798 


1.018730 9014 


0.148 x 10~ 6 


0.3 


1.040 818 2207 


1.040 8184220 


0.210 x 10~ 6 


0.4 


1.070 320 0460 


1.070 320 2889 


0.242 x 10" 6 


0.5 


1.106 530 6597 


1.106 530 9344 


0.274 x 10~ 6 


0.6 


1.148811 6361 


1.148 8119343 


0.298 x 10" 6 


0.7 


1.196 585 3034 


1.196 585 6186 


0.314 x 10~ 6 


0.8 


1.249 328 9641 


1.249 329 2897 


0.325 x 10" 6 


0.9 


1.306 569 6598 


1.306 579 9912 


0.331 x 10~ 6 


1.0 


1.367 879 4412 


1.367 879 7744 


0.333 x 10" 6 
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Methodes de Runge-Kutta Methode de Runge-Kutta d'ordre 4 



On semble avoir LA methode. 

Est-il possible qu'il soit avantageux d'utiliser une methode moins 
precise (demandant moins d'effort de calculs a chaque iterations) avec 
un pas de temps plus petit plutot qu'une methode tres precise 
(couteuse d'un point de vue calculatoire) et demandant moins 
d'iterations? 

On voudrait pourtant tenir compte de tout les facteurs : 

9 nombre d'iterations pour une valeur donnee de t Q 

9 quantite de calcul 

9 erreur sur I'approximation 
On mesure le cout du calcul par le nombre total devaluation de f 

9 RK4 demande 4 evaluations de f par iteration 

9 RK2 demande 2 evaluations de / par iteration 

9 Taylor d'ordre 2 1 evaluation de /, et de ces 2 derivees partielles 

• Euler explicite 1 evaluation de f 
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Methodes de Runge-Kutta Methode de Runge-Kutta d'ordre 4 



En utilisant y' = — y + 1 + 1 comme illustration on obtient en mesurant 
I'erreur en y(1) 



Schema 


h 


# de pas 


# eval de / 


erreur 


Euler explicite 


0.025 


40 


40 


0.464 x 10" 2 


RK2 Euler mod. 


0.05 


20 


40 


0.159 x 10" y 


RK4 


0.1 


10 


40 


0.333 x 10" b 



Conclusions : 

• II est toujours mieux d'utiliser la methode dont I'ordre est le plus 
eleve possible 

9 On comprend pourquoi RK4 est la plus utilisee. 
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La forme generale d'un systeme de m equations differentielles avec 
conditions initiales s'ecrit : 



f y\(t) = Mt,y,(t),y 2 (t),--,y m (t)) 
/ 2 {t) = W,y^t),y 2 (t),---,y m (t)) 

/s(0 = W,yi(t),Mt),'-,ym{t)) 



(yi(fc) = yi,o) 
(y2(fc) = Yw) 
Mb) = y&,o) 



I / m (0 = W,yi (0,^(0, •••.y m (0) (y m (fc) = ym,o) 

lei encore, on note yi{t n ), la valeur exacte de la /"• variable dependante 
en t = t n ety /jn , son approximation numerique. 

En general les equations sont couplees : il y a dependance des 
inconnues entre elles. 

Les m conditions initiales a droite nous garantisse I'unicite de la 
solution. 
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D'une maniere compacte, on ecrit un systeme sous la forme vectorielle 



dt 



(t) = f(t,y(t)) 



9 Toutes les methodes s'etendent presque sans changement dans 
le cas vectoriel. Comme pour la rsolution de systeme non lineaire 
on n'a qu'a faire un "passage" en dimension superieure dans les 
schemas deja construit. 
Par exemple Euler explicite devient : 

y„+i = y„ + hf(t n , y n ) 

9 La condition initiale devient vectorielle y(t Q ) = y 

9 Les considerations sur la precision s'appliquent encore : RK4 est 
a favoriser. 
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La forme generale d'une equation differentielle d'ordre m avec 
conditions initiales est : 

y^ m) (t) = f(t,y(t),y^(t),y (2) (t), . . . ,/ m -'\t)) 

ou yW(f) designe la /"■ derivee de y(t). 

Pour avoir I'unicite il nous faut m conditions initiales (la solution 
comportera m constantes que Ton determine de maniere unique avec 
m valeurs connues (et independantes) de y). Les conditions initiales 
sont donnees par 



(m-1). 



y(t ) = d, y^(k) = c 2 , • • • ,y [m -'Kk) = c m . 

On a des outils de resolution pour les systemes d'equations 
differentielles d'ordre 1 . On va transformer notre equation d'ordre m en 
m equations d'ordre 1. 
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Le principe de base consiste a ecrire I'equation d'ordre superieur sous 
la forme d'un systeme de m equations differentielles 



< 


f >i(0 

y£(0 


= y 2 (0 
= y 3 (0 
= y 4 (0 








/■■(to) 

y2(^o) 
ys(fe) 


= 0\ 

= c 2 
= c 3 




yLi(0 

I / m (0 


= y m (0 
= f(t,y: 


(t),y 2 (t),-- 


■ ,y n 


7(0) 


y m -i(fc) 

yA77(?0) 


= C m -1 


on a 


i bien 















A™) i 



A m ~), 



,0) 



yr(0 = ** (0 = ■••>% (0 = ^.yi(0,y 2 (0,--- ,y m (0) 
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Soil une equation differentielle : 

y'(t) = f(t,y(t)) 

Le principe de ces methodes consiste a integrer I'equation dans 
I'intervalle [t n , t n+ -\], c'est-a-dire : 



y'(t)dt= / f(t,y(t))dt 

tn Jtn 



ou encore : 

ftn+1 

y(Wi) = y(t„) + f(t,y(t))dt 

Jtn 

Cela nous amene a un algorithme de la forme : 



y n+ :=y n + f(t,y(t))dt 

Jtn 
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A partir d'une table de differences divisees, on remplace f 

f(t,y(t)^ Pm (t) 

par un polynome d'interpolation. 
Ainsi I'algorithme devient 

y n +\ =y n + f(t, y(0) <*«y/i+ / Pm(t) dt 

Jtn Jt„ 

Si on choisit les noeuds t n , f n _i , f n _ 2 , . . . , on obtient la famille des 
methodes d'Adams-Bashforth. 
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Formules d'Adams-Bashforth 


yn+^ -- 


= Yn + hfn (ordrel) 


Yn+\ -- 


= y n +2$ f n- f n-\) (ordre2) 


Yn+\ -- 


= y n + j2^ 23fn ~ 1 6/r "" 1 + 5/ "- 2) (ordre 3) 


Yn+\ -- 


= Yn+ 24(55/n - 59 f n _ : + 37f n _ 2 - 9f n _ 3 ) (ordre 4) 
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» La methode d'Adams d'ordre 1 correspond a la methode d'Euler 
qui est un schema a un pas. 









(RK2). 
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» La methode d'Adams d'ordre 1 correspond a la methode d'Euler 
qui est un schema a un pas. 

9 La methode d'Adams d'ordre 2 est un schema a 2 pas car depend 
de y n et de y n _i . 



t—l _.ui t 



(RK2). 
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La methode d'Adams d'ordre 1 correspond a la methode d'Euler 
qui est un schema a un pas. 

La methode d'Adams d'ordre 2 est un schema a 2 pas car depend 
de y n et de y n _i . 

Pour demarrer un tel schema, il faut connaTtre y et y-\ . On peut 
utiliser une methode a un pas du meme ordre pour calculer y-\ 
(RK2). 

En qeneral, la methode d'Adams d'ordre m est un schema a m 



Maths-Info (ENSAIV 



La methode d'Adams d'ordre 1 correspond a la methode d'Euler 
qui est un schema a un pas. 

La methode d'Adams d'ordre 2 est un schema a 2 pas car depend 
de y n et de y n _i . 

Pour demarrer un tel schema, il faut connaTtre y et y-\ . On peut 
utiliser une methode a un pas du meme ordre pour calculer y-\ 
(RK2). 

En general, la methode d'Adams d'ordre m est un schema a m 
pas. 
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Formules d'Adams-Moulton 


y n +^ -- 


= y n + hf n+1 (ordrel) 


y n +i = 


= Vn + gCWi + fn) (ordre2) 


y n +i = 


= y n + ^ (5^+1 + 8f„ - f n ^ ) (ordre 3) 


y n +\ -- 


= /" + 24 ( 9/r "+ 1 + 1 9/ " ~ 5/ "- 1 + f "-2) (° rdre 4 ) 
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» Ces methodes sont basees sur les noeuds t n+ -\ , t n , ? n -i , t n -2, ■ ■ ■ 
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» Ces methodes sont basees sur les noeuds t n+ -\ , t n , ? n -i , t n -2, ■ ■ ■ 
9 A cause de la presence de y n+1 dans le membre de droite, le 

schema devient implicite, i.e. exige la resolution d'une equation 

non lineaire pour calculer y n+1 . 



Maths-Info (ENSAIV 



9 Ces methodes sont basees sur les noeuds t n+ -\ , t n , ? n -i , tn-2, ■ ■ ■ 
9 A cause de la presence de y n+1 dans le membre de droite, le 

schema devient implicite, i.e. exige la resolution d'une equation 

non lineaire pour calculer y n+1 . 

9 Pour contourner cette difficulty, on combine les formules 
d'Adams-Bashforth et d'Adams-Moulton en des schemas dits 
predicteurs-correcteurs. 
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Schemas de prediction-correction 



Xn+1 = yn + hfn 
y n+1 = yn + htf + i 



(ordre 1) 



r m 



Yn+n^fn-fn-O 



y n+1 = y n + -{f n ' + f„) (ordre 2) 
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Schemas de prediction-correction (suite) 



y£ +1 = Yn + ^Wn-Wn^+bfn^) 



12 



y n+1 = y„ + ^ (5^ +1 + 8f„ - f„_i ) 



(ordre 3) 



y P n+: = Yn + ^r(55f„ - 59/ n _i + 37/„_ 2 - 9f„_ 3 ) 



24 



y n+1 = y„ + — (9/£ +1 + 19/,,- 5^_i + f„_ 2 ) (ordre 4) 
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Schema de prediction-correction d'ordre 2 


t 


fn 


y n 


\y(t n )-y n \ 


0.0 


— 


1.000 000 000 





0.1 


— 


1.004 837 500 


0.819 640 x 10" 7 


0.2 


1.019111875 


1.018 640 031 


0.907 218 x 10" 4 


0.3 


1.041085 901 


1.040 653 734 


0.164 486 x 10" 3 


0.4 


1.070487 675 


1.070 096 664 


0.223 381 x 10~ 3 


0.5 


1.106614851 


1.106 261088 


0.269 571 x 10" 3 


0.6 


1.148 826 758 


1.148 506 695 


0.304 940 x 10- 3 


0.7 


1.196 543 746 


1.196 254173 


0.331 129 x 10" 3 


0.8 


1.249 241382 


1.248 979 396 


0.349 568 x 10~ 3 


0.9 


1.306445195 


1.306 208166 


0.361 493 x 10" 3 


1.0 


1.367 725 911 


1.367 511462 


0.367 979 x 10" 3 
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Rappels 



Un schema numerique a un pas est une equation de recurrence de la 
forme : 

y n +i = y n + h(f)(t n ; y n ; h) 

fn+1 = t n + h 



'erreur de consistance au pas n est par definition I'erreur commise sur 

y„ + i, lorsqu'on prend pour les valeurs precedentes des y^ les valeurs 

exactes z(fo), ce qui donne la definition suivante pour les methodes a 

un pas : 

L'erreur de consistance est la suite 

e n = z(t n+: ) - z{t„) - h4>{t n ; z(t„); h) 
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